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An exact solution is found for the field of a dipole over plane, finitely conducting earth 
through an atmosphere in which the refractivity varies exponentially with height. The 
method of Hanke] transforms is used, so thai the final result, takes the form of an integral 
with infinite limits. This integral is evaluated numerically for a typical profile and for sev- 
eral wavelengths, and the results are compared with those for a homogeneous atmosphere. 
At short wavelengths and large distances, the stratified atmosphere above a plane earth can 
act as a very efficienl waveguide. 



1. Introduction 

The solutions of Maxwell's equations for a homo- 
geneous, nonconducting atmosphere over both plane 
and spherical finitely conducting earth have been 
known for many years. Almost two decades ago 
it was revealed by experiments that the field strength 
at very short wavelengths and large distances is 
greater than that calculated from "airless-earth" 
theory by tens of decibels. The favorite theory for 
explaining this excess is the scattering from randomly 
located irregularities of refractivity in the air lying 
within the beams of the two antennas [Booker and 
Gordon, 1950]. Other authors have maintained 
that the strong transhorizon fields can be explained 
only by postulating a prominent contribution from 
the ever-present stratification [Carroll and Ring, 
1955]. An exact solution for a smooth spherical 
earth and a smooth, continuous profile of refrac- 
tivity (e.g., exponential) is inordinately difficult. 
Much work has been done on the basis of solving the 
wave equation by WKB methods, [Bremmer, 
1960], and also by an extension of the latter to the 
rigorous solution [Bremmer, 1962]. Other authors 
have employed quasi-optical methods in an effort 
to allow for the effects of stratification [Bullington, 
1963]. A somewhat related problem involving 
reflection from a stratified ionosphere has been dis- 
cussed [Wait and Walters, 1963a and 1963b]. 

It is our opinion that the problem can be formula- 
ted in such a manner, especially in the case of 
"smooth" profiles, that the mathematical difficulties 
can be overcome with the aid of digital computers. 
This paper illustrates such a method in the relatively 
simple case of a plane earth and an exponential 
profile. 



2. Maxwell's Equations and Boundary 
Conditions 



A harmonic time function 



is assumed, and 



the MKS system of units is employed. Then for 
a source-free region, Maxwell's equations become 



curl HWweE 
curl E= — /co/iH 
diveE=0 
div H=0 



(1) 



Following Bremmer [1949], we assume a vertical 
electric dipole and express field strength in terms 
of a scalar Hertz potential II: 



E=(k /k 2 ) curl curl (a 2 Hl), 
K=i(ko/o)fi) curl (a 2 Ml), 



(2) 
(3) 



where a 2 is a unit vertical vector, k 2 =u 2 iJLe, and 
kl = u> 2 ix^. (Zero subscripts designate vacuum val- 
ues.) Assume m = Mo> e=e(z). Substitution of (2) 
and (3) into (1) yields the scalar wave equation which 
must be satisfied: 



where 



v 2 n+&* 2 n=o, 
~ k d?W 



2 =k' 



(4) 
(5) 
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It is readily shown that, for a smooth profile of 
tropospheric refractivity, the second term in the 
right-hand member of (5) is negligible for frequencies 
greater than about 50 kc/s. For simplicity, we 
assume k*=k. 

For an elementary electric dipole located at the 
point z=h, r—0, (cylindrical coordinates), (4) 
becomes valid everywhere if we write 



where 



v'n+Pn^s, 



S=£i(r)*(*-h), 



(6) 



.Idl 



coe 

The solutions are subject to the boundary condi- 
tions that the tangential components of E and H 
must be continuous at the boundary between air 
and earth, and that the behavior at great height 
is that characterizing free-space propagation. When 
a refractivity profile is chosen with all derivatives 
continuous, the latter requirement will be satisfied 
by proper choice of the solutions of the differential 
equation in z. The former requirement leads to the 
equations ^ 

E rl = E r2 I 

f(z=0) (7) 

where 1 refers to air and 2 to earth. The com- 
ponents of E and H in terms of II are found from 
(2) and (3). Specifically 






h b 2 (ku) ' 

k 2 drdz 
h d(M) 



oo/j, dr 



(8) 



Therefore, the boundary conditions specified by (7) 
require that 



1 &(* i n 1 )_l d(fc 2 n 2 )' 



k\ dz a 
k 1 Ui=k 2 U 2 



dz 



■(2 = 0). 



(9) 



It is necessary to find a solution of (6) which also 
satisfies (9) and the radiation condition at infinity. 
This is a unique solution. 

3. Solution of the Wave Equation 

In this paper, we assume a smooth exponential 
profile of refractivity such that 



k 2 (z)=k 2 (l+ae- z/H ), 



(10) 



i More precisely, the wave equation is 

(vH-fc* 2 )n=s 



Idl 



2ttT(jO V ' 



— 8(r)d(z-h). 

e 0« 



where a=2N s X 10~ 6 and H is the scale height. We 
first find appropriate independent solutions of the 
homogeneous equation 



(v 2 +F)n=o. 



(ii) 



The separation of variables is achieved by appli- 
cation of the Hankel transform [Sneddon, 1951] 

/» 00 

S(p, z) = I u(r 9 z)J (pr)rdr. (12) 

This leads to an ordinary differential equation in z 9 
which can be solved by any available method. Hav- 
ing found u(p, z), we can determine the original 
function by applying the inverse transform. 



u(: 



r,z)=\ u(p,z)J (pr)pdp. 



(13) 



To apply (12) to the wave equation (11), we 
multiply by rJ (pr) and integrate the result from 
to oo. The procedure requires integration by 
parts, recognition of the fact that the functions 

rJo(pr) :r- and prJ r (pr)u must vanish at r=0 and 

r=c°, and use of the Bessel equation of zero order. 
The result is, for ^=11, 



g+(F-P*)n= 



=0. 



(14) 



It may be noted that the parameter p can be loosely 
identified with k sin 0, where 6 is the angle between 
the wave normal and 2-axis. The condition p 2 ^>k 2 
corresponds to the ray-theory concept of total in- 
ternal reflection. The solutions of (14) should have 
the nature of progressive waves when p 2 <Ck 2 and of 
standing waves when p 2 ^>k 2 . It is easy to derive 
Sommerfeld 7 s solution for a homogeneous atmosphere 
from (14) and (13) by assuming k constant. 

When the exponential profile (10) is assumed, (14) 
becomes 



^+(A*+B*e-*' B )U=0, 



where 



A*=kl-p\ 
B 2 =kta. 



(15) 

(16) 
(17) 



By the change of variable v=e~ z/2H , (15) is converted 
to BessePs equation of imaginary order. We select 
the solutions 

f( P ,z)=J v tt), (18) 



ff(p,«)=r,(r), 



(19) 



where v=i2HA, ^=2HBe~ z/2H . It can be shown 
readily that only /has the proper behavior as z-^^o . 
A solution of this type was first discussed by Elias 
(1931). 

Consideration of the branch point at p=k Q leads 
to the conclusion that, on the appropriate Riemann 



1194 



sheet [Brekhovskikh, 1960], 

A branch point at p=& 2 =w&o produces no difficulties. 
Note that we define 



n 2 =(k 2 /k ) 2 =e r -im(T\ 



(20) 



where e r and o are the dielectric constant and con- 
ductivity of the earth. 

When the Hankel transform (12) is applied to the 
source function S, we obtain 



S=~5(z-h). 



(21) 



Since the boundary conditions (9) do not contain r 
explicitly, they become simply 



1 d n - N 1 d n ~ . 

MTz^^MTz^ 



(2 = 0). 



(22) 



If the earth is assumed homogeneous, (14) yields 
the solution 



where 



Ta=« m , 



y 2 =kl—p 2 =n 2 kl—p 2 . 



(23) 



Then -j^=iyu 2 , and equations (22) reduce to the 
single equation 

3gS (fclil)= §* lfil ( ^ 0) ' (24) 

We can drop the subscript on the symbol n and write 
dii 



dz 



Dii 



where 



j-^ • /ii -L ClrC\ 

/b2 A?i (tZ j 



(2 = 0). 



(25) 



Setting ki=k except in the derivative and using 
(10), we get 



w 2H 



(26) 



The second term in the right-hand member can be 
neglected except for very long wavelengths. Finally, 
it is found that the integration with respect to p 
suggested by (13) yields results only for p very close 
to k , so that we can employ the following condition 



with negligible error: 



D= 



iko^n 



(27) 



This approximation is equivalent to the introduction 



of a surface impedance, [Wait, 1962]. 

We can now combine the results given by (18), 
(19), (21), and (25) by using the Lagrange method 
of variation of parameters [Ince, 1944]: 

n(3)=^f(z)£g(3')S(z')dz'+g(z)jy(z')S(z')dz f 



Dj{0)-f 



where z r is a variable of integration and the Wron- 
skian W is given by 



W[f{z),g{z)]- 



kH 



(29) 



Evaluation of the integrals in (28) gives different 
results, depending upon whether z^>h or z<CJi. For 
z>h, 



h{p,z)= 



^/(2)[.'/W- 



D<l(0)-<)'(0) 
~Df(0)-f(0) 



./'(/')]• 



(30) 



For z<^h, z and /;, are reversed. Finally substi- 
tution in (13) yields the solution, for z>A, 



u(r, z) = 



[ g(/t) ~ g/(!)-/'(o) /(/t) ]' /o(p ' v/p ' (31) 



It is understood that / and g are functions of p as 
well as z. Since p appears in the expression for v y 
the order of the Bessel functions in (18) and (19), it 
is seen that this order ranges from large imaginary 
values to zero, then through positive real values to 
infinity as p goes from to oo along the real p-axis. 
The integration with respect to p can be thought of 
as an addition of rays propagated at various angles 
relative to the 2-axis. Small values of p correspond 
to rays which escape from the atmosphere. The 
condition p=k c corresponds to the onset of total 
internal reflection at a great height. As p continues 
to increase, the rays are reflected at lower elevations. 

It may also be noted that the integrand in (31) 
develops poles at values of p where Dj(0) =/' (0) . 
These poles lie on the positive real p axis when D is 
real or zero and slightly below this axis when D is 
complex. These poles are associated with the physi- 
cal concept of multiple reflection of rays between the 
atmosphere and the earth. For long wavelengths 
there may be but a single pole, whereas for short 
wavelengths there may be many poles. 

No serious effort has been made to evaluate the 
integral in (31) by such familiar means as the saddle- 
point method. Much better control over the accu- 
racy of the results can be achieved by numerical 
integration. Also, it is not considered very impor- 
tant at this stage to investigate the effect of antenna 
height variations. Of much more interest is the 
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variation of field strength with horizontal distance. 
Accordingly, we let z=h=0, and obtain the relatively 
simple result 



n(r,o)=-^J o "/(0) 

2xJo 



/(0)ff , (0)-/ / (0)ff(0) 



Z?/(0)-/'(0) 



J (pr)pdp 



J (pr)pdp 



(32) 



We note that 

/' (0)//(0) = - BJl (2HB)/J e (2HB) , 

p=i2H y /W Z ? P<h 



=2HJ P i -kl 
It is interesting to note that 



o>k . 



lim [/' (0)//(0)] = -B 2^5= -^f7=kT, 



p->°° 



and that the limiting form of (32) is the familiar 
Sommerfeld solution for a homogeneous atmosphere. 
Since 2Hk is normally a large number, 



Iim[/ / (0)//(0)]«-iV«-P 8 , 

and a similar coincidence with the Sommerfeld solu- 
tion results. In fact, for a homogeneous medium, 
(32) takes the form 



n( ''°>=£I' 



J (pr)pdp 
D+v/2H 



(33) 



Dimensionless forms can be obtained if welet p=k v. 
Equation (32) now becomes 



H(r, 0) 



Ok f °° 
= arj (D/fc )+V5j;(2Hi?)/J F (2H5) # 



J (k rv)vdv 



where 



(32a) 



v=i2Hk^l-v 2 v<\ 
=2Hk ^W-i v>l. 



A corresponding form is obtained for (33). Now 
only the coefficient, Ck , has dimensions. If the 
result obtained from (32) or (33) is divided by the 
free-space value 

TT 

n =— - — e 

47rr 



-ik r 



(34) 



a dimensionless ratio is obtained, which can be 
plotted as a function of r. The shape of such a 
curve is essentially the same as if the vertical com- 
ponent of E had been computed and compared with 
its free-space value. 



It is often convenient to use complex variable 
theory to evaluate the integral. It is first necessary 
to extend the path of integration along the entire 
real axis by a familiar transformation of the inte- 
grand [Sommerfeld, 1949]. Equation (32a) takes 
the form 



H(r,0>=f£ 



H^ (k rv)vdv 



(D/h) +^J' V {2HB)IJ V {2HB) 



If we divide (35) by (34), we get finally 

H {2 o> (k rv)vdv 



n f ° 

— =—k re l1c ° r 



{Dlh) +*JaJ' p (2HB)/J,(2HB) 



(35) 



(36) 



One method of evaluating the integral is to choose 
a path along the real axis. Then we use (32a), 
employ an asymptotic form for J similar to (38), 
let r be an integral number of wavelengths, and 
obtain 



n *v t Jo 



cos (k rv—Tr/4)vdv 



(Djh) +4aJ' v (2HB)IJ v (2HB) 



(37) 



Since the cosine term varies rapidly with v when 
k r becomes large, it is very helpful to choose a step 
size Vfl such that all of the integrand except the 
cosine term can be accurately represented by a 
second-degree polynomial, then integrate analyti- 
cally. In this way, the cosine term may execute 
many cycles within one step without degrading the 
accuracy of the computed value of the integral. 

A second method is based upon Cauchy's integral 
theorem. The integrand in (36) is an analytic 
function of v except at the poles, which are located 
in the second and fourth quadrants of the complex 
v plane. The path of integration is distorted into 
the lower half -plane (fig. 1), where the integrand 
vanishes sufficiently far off the real axis. However, 
the path must be taken around the branch points 
at v=l and v=n. The contribution from the 
second branch point is negligible for normal values 
of n. The cut from the branch point at v = l is 
drawn downward, and the path of integration is 
collapsed onto opposite banks of this cut. The 
total value of the integral is the sum of this branch- 
cut integral and of the sum of the residues multiplied 
by 27ri. 

Let Hi be the contribution from the branch point 
at v=l, and let #=1— ^ along this part of the path. 
The integrand on the two banks of the cut differs 
only in that v, the order of the Bessel function, 
changes sign. In the definition v=i2Hk Q ^Jl—v 2 = 
i2Hk -yjb 2j ri2bj the radical must be taken with a 
positive sign on the left bank and a negative sign 
on the right bank of the cut. For simplicity in 
writing a formula, we always take the positive sign 
and change v to -v on the right bank. Now, if we 
take advantage of the fact that the significant 
range of v is in the vicinity of 1 and decide to keep 
the minimum value of r about 10 wavelengths, we 
can use a simple asymptotic form for the Hankel 
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Figure 1. Path of integration in complex v-plane for (36). 



function: 



4 



irk r 

IT 

•k Q r 



-i(k( } rv-w/'l) 



e -i(k r— w/i) e ik Q r(l-v) ^ 



(38) 



Equation (36) leads to the form 

5l« i J— e 1 * 1 * C er k ^[ ]d8 (39) 

IIo V 7T Jo 

where the brackets are to be filled by the expression 



(P/k Q )+JaJ-' v (2HB)/J-,(2HB) 



(D/ko) + -JaJ',(2HB)/J 9 (2HB) 



This integral converges rapidly when k r is large. 

The location of a pole may be defined by the 
following relation which corresponds to a zero of 
the denominator of the integrand: 



where 



J,(2HB) = tJ'(2HB), 



T = % 



.- V! 



an' 



yjn 2 



(40) 
(41) 



Since r is normally small, the poles are located 
just below the positive real v axis, and near the 
real-axis zeros of J V {2HB). Furthermore, the real- 
axis zeros must lie in the range 0<^v<C2HB, i.e., in 
the range where v is slightly greater than 1. Equa- 
tion (40) may be solved for v s by an iterative 
method. Then 



v s =^l+v 2 J4H 2 k 2 . 



(42) 



Let n 2 be the contribution of the poles. This is 
made up of 2tt% times the residues at these poles. 
The residue at a pole is obtained by differentiating 
the denominator of (36) with respect to v, then 



substituting v=v s . We obtain 

~ (dcn()ininator) = ^^^| ; [j;(2H5)/J F (2£r5)] 
(2Hk ) 2 -y/av d 



l—v,= 



-vl 



Il2__. 

n ~ ^ (2Hh 



yiir/A 



dv 

[J v (2HB)IJ'„(2HB)], 

(2ffi ) 2 (l+^,r 2(2/^ ) 2? 
-t 2 l 2hr 

■o) 2 v«V » 

v s exp [-ihrp 2 8 /2(2Hk y} 
s=l ~ [J v (2HB)/J' v (2HB)) v=Cs (43) 

Since v s normally has only a small negative imagi- 
nary component, the damping factor for the residues 
lends to remain near unity, even for large values of 
& f. This means small attenuation of the trapped 
modes. 

The total Hertz potential relative to its free-space 
value is the sum of (39) and (43). For long wave- 
lengths, where 2TIB is small, there may be only one 
root of eq (40), therefore only one pole. For short 
wavelengths, there are many roots of (40) and many 
poles. This behavior is somewhat similar to that of 
a metallic wave-guide, in which many modes can 
propagate at sufficiently short wavelengths. 

4. Computations 

Only a few comments are required with regard to 
the computation of numerical results. The evalua- 
tion of the integral in (37) or (39) by an application 
of Simpson's method is routine except for the rapid 
oscillation of the cosine function in (37). This is a 
common problem in Fourier analysis. It is generally 
somewhat easier to use complex-plane integration, 
especially when 2HB is large and J V (2HB) has many 
zeros on the real z;-axis. 

Perhaps no other functions have been so exten- 
sively investigated as the solutions of BessePs 
equation. The work of Langer [1931] and Watson 
[1944] is especially noteworthy. Since both the 
argument and order of the Bessel function are often 
large in this application, and since the order may be 
complex, asymptotic formulas must be used and care 
must be exercised in choosing the formula which is 
appropriate to the particular region of the complex 
plane (Stokes phenomenon). It is also to be noted 
that computation tends to be difficult in the region 
where the order and argument of the Bessel function 
become equal, i.e., the region where the coefficient of 
n in (15) becomes zero. In the present problem, 
this difficulty occurs only when integration along the 
real axis is employed, since the argument of the 
Bessel function is always real. It can be avoided in 
the simple case of zero antenna heights, because a 
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familiar recursion formula can be used for computing 
J'jJ. A useful form is the following: 






H-i | «/>+i(f) 

f ~ l V, + i(r) 



(44) 



The ratio J' '/J is first computed for a larger order 
v-\-n, where n is an integer, from a simple asymptotic 
formula, and an iterative procedure is used to reach 
the final result. The location of the poles v s in 
accordance with (40) also presents some difficulties. 
The general procedure is to obtain a first estimate of 
v s and then make corrections by an iterative method 
until (40) is satisfied. Two cases may be distin- 
guished, viz, 2HB small and 2HB large. In the 
former case a first estimate of v x can be obtained from 
the formula 



J'v&HB) 



HB 



r J Vl {2HB) 2HB Vl +1 



(45) 



In the latter case, numerous poles exist, so that 
\vi\<C<C^HB, and we can use a simple asymptotic 
formula for a first estimate: 



= J^(2HBY 



,-tan(2HB-^-|) 



(46) 



The iterative method for improving upon such an 
estimate is based upon the numerical determination 
of the derivative of JjJ' with respect to v. In this 
process, the derivative appearing in (43) is also 
computed. Where several poles occur, extrapola- 
tion methods can be used to good advantage for a 
first estimate after the pole of lowest order is located. 
In the following numerical examples, both real- 
axis and complex-plane integration have been em- 
ployed, and the results were found to agree. 

5. Numerical Examples 

(a) /=100kc/s,X =3000m,e r =10,o-=10- 2 mho/m, 
^ 2 =10-il800, q= 6X10- 4 , #=6000 m, 2HK =8t= 
25.1327, 2#5=87rV6X10- 4 -:0.615624. 

Vertical electric dipole on the surface: 

t=1. 03924 exp ^0.78846, *>=i25. 1327-^1 -^ 2 . 

Location of pole: 

first estimate from (45) . . . ^ — 0.535—^0.390, cor- 
rected location . . . ^i=0. 53581— 10. 38902. 

The ratio of the Hertz potential to that in free 
space was computed with the aid of a large digital 
computer as a function of horizontal distance from 
the source. The results are shown in figure 2. The 
effect of the stratified atmosphere is discernible at 
long distances but is too small to be important 
practically. 
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Figure 2. Hertz potential relative to free-space value. 
A = 3000 m. 



(b) /=10 Mc/s, A o = 30 m, e r =10, o-=10" 2 mho/m, 
n 2 =10— il8. a=6X10~ 4 , #=6000 m, 2HK 
= 800tt = 2513.27, 2HB=800wyfoX 10- 4 =61.5624. 

Vertical electric dipole on the surface: 



r = 0.11243 exp i 1.06065, v=i 2513.27 Vl 



— Tf 



Location of poles: 

first estimate from (46)— ^» 1.6565— i 0.06247, 

corrected location— v x = 1.6697— i 0.063495. 



v 2 = 3.7271- 
^ 3 = 5.8313- 
p 4 = 7.9857- 
^5=10.1942- 
^ 6 =12.4616- 
* 7 = 14.7932- 
» 8 =17.1954- 
v 9 = 19.6762- 
^10 = 22.2447- 



i 0.064808 
i 0.066138 
* 0.067486 
i 0.068856 
i 0.070251 
i 0.071672 
i 0.073123 
i 0.074610 
i 0.076132 



^ii = 24.9130- 

*>i 2 = 27.6957- 
^13 = 30.6126- 
^14 = 33.6901- 
^15 = 36.9658- 
*/ 16 = 40.4966- 
^ = 44.3770- 
^i 8 = 48.7895- 
z,i 9 = 54.2070- 



-i 0.077705 
■i 0.079327 
■i 0.081015 
■i 0.082779 
■i 0.084639 
■i 0.086638 
■i 0.088816 
■i 0.091324 
■i 0.094847 



The results in the form of relative Hertz potential 
versus distance are shown in figure 3. In this case, 
the effect of the stratified atmosphere is very prom- 
inent at long distances. It is apparent that the 
trapped modes account for a slowly decaying field 
in this distance range. The fluctuations due to 
mode interference become so numerous and irregular 
for distances greater than 1500 km that only the 
square root of the sum of the squares of mode 
amplitudes was plotted. In this region, the ampli- 
tude is a very sensitive function of distance or of 
variations in the parameters of the refractivity 
profile. Temporal variations of these parameters 
can be expected to produce fading at a fixed distance. 
The situation is somewhat similar to that believed 
to exist in short-wave transhorizon propagation. 
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Figure 3. Hertz potential relative t o free-space value. 
A„=30m. 

6. Conclusions 

A practical method has been developed which 
allows computation of the Hertz potential produced 
by a vertical dipole over a finitely conducting earth 
and in an atmosphere with an exponential profile of 
refractivity. The mathematical solution is an in- 
tegral, which must be evaluated by numerical 
methods. The Caucliy integral theorem may be 
applied to show the existence of one or more trapped 
modes. For long wavelengths, such as 3000 m, the 
attenuation of the one trapped mode is rather large, 
and the influence of the stratification of the atmos- 
phere is slight. In the range of 10,000 km and 
beyond in figure 2, the result comes almost entirely 
from the branch-cut integral (39). For shorter 
wavelengths, such as 30 m, numerous trapped modes 
exist and the attenuation as a function of distance is 
amazingly low. Beyond 600 or 700 km in figure 3, 
the field results in large part from the residues (43) . 
At long distances, these residues interfere to produce 
rapid fluctuations of field strength as a function of 
distance. This could lead to fading at a fixed dis- 
tance if the atmospheric structure varied slightly 
with time. The situation seems somewhat analogous 
to that encountered in scatter propagation. 

This paper deals with a case which is not in keep- 
ing with the spherical geometry of the real world. 



However, the wave equation can be solved readily 
for this case, whereas the same is not true in a 
spherical system. It is believed that certain of the 
concepts developed here may be useful in an attack 
on the more difficult problem of propagation in a 
stratified medium surrounding a sphere. 
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